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FINAL  REPORT: 

DENSE-GAS  DISPERSION  IN  COMPLEX  TERRAIN 

D.L.  Ermak 

Lawrence  Livermore  National  Laboratory 
Livermore,  California  94550 

ABSTRACT 

A  dense-gas  version  of  the  ADPIC  Lagrangian  particle,  advection-diffusion  model  has  been 
developed  to  simulate  the  atmospheric  dispersion  of  denser-than-air  releases  over  complex 
terrain.  In  developing  the  model,  it  was  assumed  that  the  dense-gas  effects  could  be  described  in 
terms  of  the  vertically-averaged  thermodynamic  properties  and  the  local  height  of  the  cloud.  The 
dense-gas  effects  were  treated  as  a  perturbation  to  the  ambient  thermodynamic  properties 
(density  and  temperature),  ground  level  heat  flux,  turbulence  level  (diffusivity),  and  windfield 
(gravity  flow)  within  the  local  region  of  the  dense-gas  cloud.  These  perturbations  were  calculated 
from  conservation  of  energy  and  conservation  of  momentum  principles  along  with  the  ideal  gas 
law  equation  of  state  for  a  mixture  of  gases.  ADPIC,  which  is  generally  run  in  conjunction  with  a 
mass-conserving  wind  flow  model  to  provide  the  advection  field,  contains  all  the  dense-gas 
modifications  within  it.  This  feature  provides  the  versatility  of  coupling  the  new  dense-gas 
ADPIC  with  alternative  wind  flow  models.  The  new  dense-gas  ADPIC  has  been  used  to  simulate 
the  atmospheric  dispersion  of  ground-level,  colder-than-ambient,  denser-than-air  releases  and  has 
compared  favorably  with  the  results  of  field-scale  experiments. 
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A  popular  approach  to  simulating  the  atmospheric  dispersion  of  trace-gas  releases  under 
realistic  conditions  of  terrain  and  varying  winds  is  to  use  a  Lagrangian  particle,  advection- 
diffusion  model  such  as  ADPIC1- 2  These  codes  are  well  suited  to  operational  uses  for 
emergency  response  and  accident  preparedness  purposes3  due  to  their  robustness  and  relatively 
high  computational  speed  on  present  day  computers.  While  advection-diffusion  models  have 
been  extensively  used  and  validated  against  numerous  field  tests  of  trace  gas  releases,  they  do  not 
contain  the  necessary  physics  to  simulate  the  atmospheric  dispersion  of  denser-than-air-releases. 
To  overcome  this  limitation,  we  have  developed  a  dense-gas  version  of  ADPIC. 

The  atmospheric  dispersion  of  denser-than-air  gases  requires  a  somewhat  different  approach 
to  modeling  than  the  more  conventional  one  used  for  the  dispersion  of  trace  gases4’  In  a  trace- 
gas  release,  the  quantity  of  material  released  to  the  atmosphere  is  too  small  to  have  any 
significant  effect  upon  the  atmospheric  flow  into  which  it  is  mixing.  Consequently,  trace-gas 
dispersion  is  controlled  solely  by  the  advective  and  diffusive  properties  of  the  ambient 
atmosphere.  On  the  other  hand,  a  dense-gas  release  behaves  more  like  an  independent, 
continuous  cloud  whose  physical  properties  (density,  temperature,  turbulence  level)  are 
significantly  different  than  those  of  the  ambient  atmosphere.  Under  these  conditions,  the 
dispersion  of  a  dense-gas  cloud  is  controlled  as  much,  if  not  more,  by  these  in-cloud  properties  as 
by  the  conditions  existing  in  the  ambient  atmosphere. 

Several  major  effects  are  observed  in  the  dispersion  of  dense-gas  clouds  released  at  ground 
level  that  are  not  observed  in  the  dispersion  of  trace  emissions.  One  is  a  reduction  of  turbulent 
mixing  in  the  vertical  direction  of  the  vapor  cloud  with  the  ambient  atmosphere  due  to  stable 
density  stratification  of  the  cloud  relative  to  the  above-lying  ambient  air.  Another  is  the 
generation  of  gravity-spreading  flow  due  to  density  gradients  in  the  horizontal  direction.  These 
two  effects  result  in  a  lower  and  significantly  wider  cloud  than  is  observed  when  a  trace  or 
neutral  density  gas  is  released.  In  addition,  dense-gas  clouds  tend  to  follow  the  downhill  slope 
independent  of  the  wind  direction  and  can  become  “trapped”  in  valleys  or  low  spots.  And  when 
the  release  is  a  cold,  dense-gas  (i.e.,  LNG),  the  ground  surface  heat  flux  into  the  cloud  may  also 
play  a  significant  role  in  the  dispersion  of  the  cloud. 

Dense-gas  dispersion  effects  are  most  pronounced  when  the  ambient  wind  speed  is  low  and 
the  atmospheric  stability  conditions  are  stable.  As  the  cloud  mixes  with  the  surrounding  ambient 
atmosphere,  the  cloud  becomes  more  dilute,  the  in-cloud  properties  approach  ambient  levels,  and 
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the  above  mentioned  effects  begin  to  play  a  less  significant  role.  Eventually,  after  a  considerable 
amount  of  dilution,  the  originally  dense-gas  cloud  begins  to  disperse  like  a  trace-gas  cloud  where 
dispersion  is  primarily  controlled  by  the  ambient  wind  speed  and  atmospheric  stability. 

In  the  development  of  the  dense-gas  version  of  ADPIC,  we  assumed  that  the  dense-gas 
effects  could  be  described  by  the  vertically-averaged  thermodynamic  properties  of  the  cloud  and 
the  local  cloud  height.  This  simplification  has  been  successfully  employed  in  one-dimensional 
dense-gas  dispersion  models6  and  basically  assumes  that  the  dense-gas  cloud  is  relatively  thin  in 
comparison  to  its  horizontal  dimensions.  We  then  treated  the  dense-gas  effects  as  perturbations  to 
the  ambient  conditions  within  the  local  region  of  the  dense-gas  cloud7.  Submodels  for  the 
friction  velocity,  turbulent  diffusivity,  and  ground  heating  rate  were  adapted  from  the  SLAB6  and 
FEM3^  models,  both  of  which  have  been  tested  and  validated  extensively  by  comparison  with 
laboratory  and  field-scale  dispersion  tests9-11. 

In  this  report,  we  describe  the  dense-gas  version  of  ADPIC,  the  results  from  simulations  of 
two  field-scale  dense-gas  dispersion  experiments,  and  a  calculation  of  a  dense-gas  flow  up  a 
ramp.  The  following  section  includes  a  mathematical  description  of  the  Lagrangian  particle 
advection-diffusion  model  and  the  main  dense-gas  modifications  to  the  model.  This  is  followed 
by  a  description  of  the  code  structure  with  regard  to  the  new  modifications  and  then  the  results 
from  the  model  simulations.  The  experiments  being  simulated  involved  the  large-scale  release  of 
liquefied  natural  gas  (LNG)  onto  a  pond  of  water  with  the  atmospheric  dispersion  occurring 
principally  over  land  with  nearly  flat  terrain.  The  ramp  calculation  involves  the  dispersion  of  an 
instantaneous  sulfur  dioxide  (SO2)  release  at  the  base  of  a  ramp  and  demonstrates  the  ability  of 
the  model  to  simulate  dispersion  over  variable  terrain. 
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MODEL  DESCRIPTION 

LAGRANGIAN  PARTICLE.  ADVECTION-DIFFUSION  MODEL 

The  ADPIC  atmospheric  dispersion  model  is  based  upon  the  single  conservation  of  species 
principal.  This  principal  is  expressed  by  the  advection-diffusion  equation  as 

If*  _ 

^■+V.(t/.C)  =  V-(ffVC)  (1) 

where  C  is  the  concentration  of  the  gas  released  into  the  atmosphere,  Ua  is  the  nondivergent 
(V  •  Ua  =  0)  advection  field  (ambient  wind  velocity),  and  K  is  the  ambient  atmospheric  turbulent 
diffusivity.  The  ambient  advection  field  Ua  is  usually  calculated  in  a  separate,  mass-conserving, 
wind  flow  model  from  a  limited  set  of  wind  speed  and  direction  observations.  The  turbulent 
diffusivity  K  is  a  function  of  the  ambient  atmospheric  conditions  and  is  calculated  within  the 
ADPIC  code. 

In  the  Lagrangian  particle  numerical  integration  technique  used  in  ADPIC,  the  mass  of  the 
released  gas  is  represented  by  marker  particles.  The  dispersion  of  the  released  gas  is  simulated  by 
the  marker  particle  trajectories  which  are  calculated  from  a  displacement  equation  that  is  derived 
from  the  advection-diffusion  equation.  Consequently,  the  distribution  of  marker  particles  in  time 
and  space  represents  the  concentration  distribution  of  the  released  gas.  In  the  ADPIC  code,  the 
concentration  within  a  particular'  volume  element  is  calculated  by  summing  the  total  mass  within 
the  volume  element  as  represented  by  the  marker  particles  and  then  dividing  this  sum  by  the 
volume  of  the  volume  element 

The  marker  particle  displacement  equations  used  in  ADPIC  are  stochastic  in  nature  and 
include  both  deterministic  and  random  terms.  The  form  we  use  is  adapted  from  Boughton,  et 
a l .12  and  is 

X(t  +  At)  =  X(t)  +  Ua-At  +  Rx  {At) 

Y(t  +  At)  =  Y{t)  +  Vtt  ■  At  +  Ry (At)  (2a) 

Z(t  +  At)  =  Z(t)  +  Wa-At+vo-At  +  R.  (A/) 
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where  Ua,Va,  and  Wa  are  the  three  components  of  the  ambient  wind  velocity  and 

RX,R  ,  and  Rz  are  three  independent  Gaussian  random  displacements  with  the  mean  and  mean 

square  properties 


<*,)  =  <>,  (Rl) 

W=°.(^2> 

W-0.  (*,’} 

K 

V0 

0 

In  Eq.  (2),  the  terms  involving  Ua  are  the  advective  displacement  and  those  involving  K  are  the 
diffusive  displacement.  The  advective  displacement  is  deterministic,  while  the  diffusive 
displacement  includes  both  deterministic  and  random  terms.  When  the  diffusivity  is 
homogeneous,  as  is  assumed  here  for  the  two  horizontal  directions,  the  diffusive  displacement 
contains  only  random  terms.  However,  when  it  is  inhomogeneous,  as  is  the  case  in  the  vertical 
direction,  then  there  is  an  additional  deterministic  term  involving  the  gradient  of  the  diffusivity. 

Calculating  the  marker  particle  trajectories  using  Eq.  (2)  and  then  calculating  the 
concentration  from  the  distribution  of  marker  particles  is  equivalent  to  calculating  the 
concenU'ation  by  directly  solving  the  advection-diffusion,  partial  differential  equation  as 
expressed  in  Eq.  (1).  A  major  advantage  of  the  stochastic,  Lagrangian  particle  method  is  that  it 
does  not  require  the  use  of  a  fixed  grid  to  calculate  the  particle  trajectories,  while  the  numerical 
solution  of  a  partial  differential  equation  does.  This  advantage  is  especially  important  when 
simulating  dense-gas  dispersion  since  the  low-lying  nature  of  dense-gas  clouds  may  preclude 
(from  a  practical  standpoint)  the  accurate  determination  of  the  cloud  height  and  vertical  cloud 
structure,  when  required  to  use  a  fixed  grid. 

DENSE  GAS  EFFECTS  AND  MODIFICATIONS 

Our  treatment  of  dense-gas  dispersion  within  the  ADPIC  advection-diffusion  model 
framework  requires  modification  of  the  ambient  wind  field  U„  and  the  vertical  diffusivity  K., 
and  inclusion  of  surface  heating  Jg  of  the  cloud  at  the  ground.  The  modifications  to  the  wind 


=  2K-At 


2- Ky-  At 


1  2  •  Kg  •  At  +  Vg  ’  At‘  ,  where 

Z„  =  Z(o). 

3K,(Z.) 

-ST' and 

average  value. 
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field  and  diffusivity  within  the  local  region  of  the  dense-gas  cloud  are  due,  for  the  most  part,  to 
the  higher  density  of  the  cloud  relative  to  the  ambient  atmosphere,  and  they  directly  effect  the 
advection  and  diffusion  terms  in  Eq.  (2)  for  the  marker  particle  displacement.  Surface  heating  of 
the  cloud  is  coupled  to  these  advection  and  diffusion  effects  since  it  warms  the  cloud  and  thereby 
lowers  the  cloud  density.  Surface  heating  also  adds  a  term  to  the  vertical  displacement  equation 
due  the  expansion  of  the  cloud  associated  with  the  heat  gain. 

Veitical  Averaging.  The  dense-gas  modifications  are  calculated  by  using  a  two-dimensional 
extension  of  the  one-dimensional  SLAB6  model  approach.  This  approach  is  based  on  the 
assumption  that  the  dense-gas  cloud  height  is  relatively  low  in  comparison  to  its  horizontal 
dimensions.  Consequently,  cloud  behavior  is  assumed  to  be  described  by  the  layer-(height-) 
averaged  cloud  properties  and  the  thickness  (cloud  height)  of  the  layer. 

Using  density  as  an  example,  the  vertically  averaged  cloud  density  is  defined  to  be 

p{x,y)  =  |  •  f'dz-p{x,y,z),  (3) 

where  p  =  p{x,y )  is  density,  h  =  h(x,y)  is  the  cloud  height,  and  both  p  and  h  are  functions  of 
the  two  horizontal  coordinates  x  and  y.  For  a  cloud  at  ground  level,  the  cloud  height  h  is  defined 
in  terms  of  the  mean  squared  marker  particle  height  such  that 

h2  =  3-Z\  (4) 

For  the  puipose  of  calculating  the  dense  gas  modifications,  a  three-dimensional  description  of 
the  dense  gas  cloud  is  obtained  by  assuming  vertical  profiles  for  the  cloud  properties.  Using  the 
above  example  of  density. 


p(x,y,z)=p{x,y)-f(z),  (5) 

where  f(z)  is  a  convenient  function  of  height  such  as  a  step,  Gaussian,  or  other  simple  analytic 
function.  (The  three-dimensional  concentration  distribution  is  still  determined  from  the  marker 
particle  locations;  however,  the  concentration  distribution  is  not  required  at  each  step  of  the 
numerical  integration  of  the  marker  particle  trajectories.)  The  vertically  averaged  cloud 
properties  of  primary  importance  are  cloud  height,  density,  and  temperature.  The  various  dense 
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gas  modifications  are  all  expressed  in  terms  of  these  layer-averaged  quantities  and,  when 
necessary,  an  assumed  vertical  profile. 

Gravity  Flow.  Dense-Gas  Diffusivitv.  and  Ground  Heating.  The  dense-gas  modification  to 
the  advection  field  replaces  the  ambient  wind  field  Ua  with  a  new  wind  field  U  that  includes  the 
effects  of  gravity  flow.  Treating  gravity  flow  as  a  perturbation  to  the  ambient  flow,  the  new  wind 
field  can  be  expressed  as 

U(x,  y,z)  =  Ua(x,y,z)+  Ug  (x,y)  •  g,(z)  (6) 

where 

Ua(x,y,z)  =  ambient  wind  field 
Ug(x,y)  =  vertically  averaged  gravity  flow  field  perturbation 

gi(z)  =  velocity  profile  function  where  i  =  h  (horizontal)  or  /  =  v  (vertical). 

The  vertically  averaged,  horizontal  components  of  the  gravity  flow  field  are  obtained  by  equaling 
the  gravitational  force  terms  (using  the  hydrostatic  approximation)  to  the  dissipative  terms  due  to 
air  entrainment  and  surface  friction,  yielding 

P0)eUg  =  ~fy°'5S  ^ P ~P^~g  (P-P«)h^-Puf  U* 

.  (7) 

pO)eVg  =-— [0.5g  (p-pa)h2]-g{p-pa)h  —  -pufVg 

where  coe  =  entrainment  rate 

g  =  acceleration  of  gravity 
h  =  cloud  height  above  terrain 
H  =  terrain  height 
uf  =  surface  friction  coefficient 

The  entrainment  rate  0)e  is  the  rate  at  which  air  mixes  into  the  dense-gas  cloud  from  above  and  is 
taken  to  be  proportional  to  the  vertical  diffusivity  Kz(h)  divided  by  the  cloud  height  h.  The 
vertical  component  of  the  gravity  flow  field  is  obtained  from  the  nondivergence  condition 
(V  •  U  =  o)  and  is 
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r^+5v!'i 

dx  dy 


■ildzgh(z) 


(8) 


Equation  (8)  defines  the  relationship  between  the  horizontal  velocity  profile  function  gh(z)  and 
the  vertical  velocity  profile  function  gy(z).  By  choosing  gh(z)  judiciously,  both  g/((z)  and  g,.(z) 
can  be  analytical  functions. 

The  new  vertical  diffusivity  is  adapted  from  the  Af-theory  diffusivity  used  in  the  FEM3^ 
three-dimensional,  conservation  equation  model  for  dense  gas  dispersion.  For  a  cloud  within  the 
surface  layer,  the  new  diffusivity  is 

z  d> 

where  O  =1  +  5- z/  L  ^ 

L-'  =L;‘+2S  ■*=  (p-pj/(p-«r) 

«.  =  friction  velocity. 


This  form  of  the  diffusivity  is  obtained  from  the  well  known  similarity  arguments  of  Monin  and 
Obukhov  and  is  such  that  it  approaches  the  ambient  diffusivity  as  the  cloud  density  approaches 
the  ambient  level. 


To  account  for  ground  heating  of  a  cold,  dense-gas  cloud,  we  use  a  simple  heat-transfer 
model 


-vft-r)  00) 

where  Jg  is  the  ground  flux  of  heat  into  the  cloud,  Vfr  is  a  representative  energy-transfer  velocity, 
Tg  is  the  ground  temperature,  and  p,Cp  and  T  are  the  local  layer- averaged  cloud  density,  specific 

heat,  and  temperature,  respectively.  In  the  simulations  of  a  cold,  dense-gas  release,  we  have 
found  it  adequate  to  replace  the  ground  temperature  Tg  with  the  ambient  air  temperature  Ta. 

Thermodynamic  Properties.  The  vertically  averaged  cloud  density  is  obtained  from  the  ideal 
gas  law  and  is  calculated  from  the  location  and  mass  of  the  marker  particles.  For  a  two  species 
system  of  released  gas  and  air  at  constant  pressure,  the  ideal  gas  law  can  be  expressed  as 
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where  p,  pe,  and  T  are  the  vertically  averaged  cloud  density,  released  gas  density,  and  cloud 
temperature,  pa  is  the  ambient  air  density  at  the  ambient  temperature  Ta,  and  Ma  and  Ms  are  the 
molecular  weights  of  air  and  the  released  gas. 

The  vertically  averaged  released  gas  density  is  defined  to  be 

P.={z°>.q,)HA-h)  (12) 

where  qn  =  q  is  the  constant  mass  associated  with  each  marker  particle,  CDn  is  a  position 
dependent  weighting  factor,  h  is  the  cloud  height,  and  the  sum  is  over  all  marker  particles  within 
the  rectangular  column  of  horizontal  area  A.  The  weighting  factor  is  calculated  by  treating  each 
marker  particles  as  if  it  were  a  “square”  with  an  area  equal  to  the  grid  cell  area  and  centered  at 
the  location  of  the  marker  particle.  This  weighted-average  approach  creates  smoother  distribution 
functions  since  each  marker  particle  contributes  to  the  average  in  surrounding  cells  as  well  as  the 
one  where  it  is  centered. 

The  cloud  temperature  is  obtained  from  the  conservation  of  energy  principal.  As  shown  in 
Ermak  and  Lange13,  the  vertically  averaged  cloud  temperature  can  be  calculated  from  the  initial 
source  temperature  of  the  released  gas,  the  ambient  temperature,  and  the  location  and  “energy 
deficit”  associated  with  each  marker  particle.  The  result  is 

T  =  [(p  -P.)  ■  C„  ■  T.  +p,  ■  C„  ■  r„]  /  (p  •  C,)  (13a) 


where  the  effective  source  temperature  is 


r  =r  +(t  -T  )•€ 

se  a  \  s  a  / 


(13b) 


and  the  vertically  averaged  energy  deficit  is 


(13c) 
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The  time- varying  marker  particle  energy  deficit  e„  (t)  accounts  for  the  addition  of  heat  to  the 
cloud  at  the  ground  surface.  This  effect  is  calculated  from  a  rate  equation  for  the  marker  particle 
energy  deficit  that  is  analogous  to  the  marker  particle  displacement  equations.  The  rate  equation 
to  first  order  in  the  time  step  is13 


E„(«  +  At)«E.(0/[l  +  (V, /*)•*«]  (14) 

where  Vh  is  the  ground  surface  energy-transfer  velocity  defined  in  Eq.  (10)  and  h  is  the  cloud 
height.  The  initial  value  of  e„  ( o )  is  1.0  so  that  the  initial  effective  source  temperature  is  equal  to 

the  actual  source  temperature. 

Surface  Heating  Expansion  Displacement.  Surface  heating  of  the  cloud  not  only  warms  the 
cloud,  but  also  causes  it  to  expand.  Calculation  of  this  expansion  is  straight  forward  if  we  assume 
that  the  expansion  due  to  surface  heating  takes  place  solely  in  the  vertical  direction.  This 
assumption  is  quite  realistic  especially  for  cold,  dense  clouds  which  are  typically  much  wider 
than  they  are  high.  Using  this  assumption,  the  surface  heating  expansion  displacement  is  found  to 
be13 


A Z„  = 


[l  +  (V„/A)A»] 


(15) 


where  Z0  =  Z(t)  is  the  height  of  the  marker  particle  at  the  beginning  of  the  time  step. 

The  total  marker  particle  displacement  in  the  vertical  direction  is  therefore 

Z(t  +  At)  -  Z(t)  =  AZh  +  AZd  (16) 

where  A ZD  is  the  dispersive  displacement  due  to  advection  and  diffusion  given  by  Eq.  (2). 
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CODE  STRUCTURE 


All  of  the  dense-gas  modifications  are  included  in  one  of  twelve  subroutines  within  the 
ADPIC  code,  as  shown  in  Fig.  (1).  Of  the  twelve  subroutines,  five  are  used  to  initialize  the 
problem,  five  are  used  to  calculate  the  various  dense-gas  effects  including  the  marker  particle 
energy  deficit  and  trajectories,  one  is  used  to  set  the  time  step,  and  one  is  used  to  generate  the 
output  files.  To  make  these  changes,  we  modified  six  existing  subroutines  and  added  six  entirely 
new  subroutines. 


MAIN 


Figure  1.  That  part  of  the  ADPIC  subroutine  tree  which  includes  the  dense-gas  modifications.  A 
dashed  underline  indicates  a  pre-existing  subroutine  that  was  modified  and  a  solid  underline 
indicates  an  entirely  new  subroutine. 

All  of  the  dense-gas  subroutines  are  called  either  directly  from  the  MAIN  subroutine  or  from 
a  subroutine  one  layer  below  MAIN.  The  initialization  subroutines  are  all  called  from  INITIAL 
and  the  output  subroutine  is  called  from  WRDPOT.  The  initialization  and  output  subroutines 
perform  the  following  functions. 

INl'iXAL-calls  the  various  initialization  subroutines. 


n 


INITl-sets  the  constants  used  in  the  dense-gas  algorithms. 


INIT2-sets  several  constants  associated  with  the  thermodynamic  properties  of  the  atmosphere 

and  checks  a  number  of  the  input  parameters  to  ensure  that  they  are  within  realistic  limits. 

GRIDSET-calculates  the  physical  properties  of  the  dense-gas  source  including  the  vertical  source 
velocity. 

TSLOPE-calculates  the  slope  of  the  terrain. 

PRBINS-generates  the  output  files  of  the  vertically  average  cloud  properties  including:  density, 
temperature,  cloud  height,  mass  concentration,  mass  fraction,  and  volume  fraction. 

The  Lagrangian  particle  simulation  of  the  dispersing  dense-gas  cloud  is  performed  from  three 
subroutines  that  are  called  directly  from  the  MAIN  subroutine.  The  subroutine  VERTAVED 
deals  with  the  vertically  averaged  cloud  properties,  MOVPARTD  deals  with  the  individual 
marker  particle  properties  and  displacement,  and  TIMSTEP  sets  the  time  step  to  be  used  in  the 
marker  particle  displacement  and  energy  deficit  equations  during  the  next  time  cycle.  These 
subroutines  perform  the  following  functions  and  call  the  following  auxiliary  subroutines. 

VERTAVED-calculates  the  local  cloud  height,  a  surface  heating  rate  parameter,  and  the 
vertically  averaged  cloud  properties  including:  energy  deficit,  mass  concentration, 
density,  temperature,  and  mass  fraction.  Calls  the  SMOOTH  subroutine. 

-SMOOTH-spatially  smoothes  the  cloud  height,  surface  heating  rate  parameter,  and  vertically 

averaged  cloud  density,  temperature  and  mass  fraction.  Calculates  the  vertically  averaged 
volume  fraction  from  the  smoothed  mass  fraction.  (Spatial  smoothing  is  performed  to 
remove  spurious  gravity  waves.) 

MOVPARTD-calculates  the  individual  marker  particle  displacements  and  energy  deficits.  Also 
calculates  the  advection  velocity  and  diffusivity  at  the  location  of  each  marker  particle. 

To  assist  in  this  calculation,  MOVPARTD  calls  the  DISPERSE  and  GRAVFLOW 
subroutines. 
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-DISPERSE-calculates  the  effective  friction  velocity,  entrainment  coefficient,  and  vertical 

diffusion  coefficient  at  the  top  of  the  cloud.  Also  calculates  additional  surface  heating 
parameters. 

-GRAVFLOW-calculates  the  three  components  of  the  gravity  flow  velocity. 

TIMSTEP-overrides  the  normal  ADPIC  time  step  calculation  at  the  end  of  each  time  cycle  and 
sets  the  time  step  to  the  constant  value  specified  in  the  input  data. 
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VALIDATION  CALCULATIONS 


We  used  the  dense-gas  version  of  ADPIC  to  predict  the  vapor  dispersion  from  two  liquefied 
natural  gas  (LNG)  releases  at  ground-level  onto  water.  The  calculations  simulate  the  vapor 
dispersion  observed  in  two  field-scale  experiments,  Burro  8  and  Burro  9,  conducted  by  the 
Lawrence  Livermore  National  Laboratory  in  1980  at  China  Lake,  California.14  The  composition 
of  LNG  is  mainly  methane  which  has  a  molecular  weight  half  that  of  air.  However,  LNG  is 
formed  by  cooling  the  gas  to  1 12°K  so  that  the  source  temperature  is  approximately  one-third 
that  of  ambient  air  (~300°K).  The  net  effect  is  that  LNG  vapor  is  about  65%  denser  than  ambient 
air  as  well  as  being  much  cooler  than  the  ambient  air  temperature.  Consequently,  these 
simulations  utilized  all  of  the  dense-gas  algorithms  implemented  into  ADPIC. 

In  the  Burro  8  and  Burro  9  experiments,  approximately  25  m3  of  LNG  (~6000m3  of  LNG 
vapor)  was  released  at  a  nearly  constant  rate  over  about  one  and  one-half  minutes.  Although  the 
release  of  LNG  was  onto  water  (in  the  form  of  a  small  pond),  most  of  the  dispersion  occurred 
over  land.  The  source  characteristics  and  meteorological  conditions  for  the  two  experiments  are 
given  in  Table  1.  The  main  difference  between  these  two  experiments  was  the  meteorological 
conditions.  The  Burro  8  test  occurred  under  lower  wind  speed  and  more  stable  ambient 
atmospheric  conditions  than  did  the  Burro  9  test 

Table  1.  Summary  of  the  Burro  8  and  Burro  9  test  conditions. 


Burro  8 

Burro  9 

LNG  volume  spilled  (m3) 

28.4 

24.2 

Mean  LNG  spill  rate  (m^/min) 

16.0 

18.4 

Mean  wind  speed  at  2  m  height  (m/s) 

1.8 

5.7 

Mean  temperature  at  2  m  height  (°K) 

306. 

308. 

Atmospheric  stability  class 

E 

D 

The  ADPIC  simulations  were  conducted  using  5000  marker  particles  within  a 
1000  x  600  x  28  m  grid.  The  horizontal  grid  spacing  used  in  the  calculation  of  the  wind  field,  the 
three-dimensional  concentration  disU'ibution,  and  the  vertically  average  cloud  properties  was 
Ax  =  25  m  and  Ay  =  15  m.  The  vertical  grid  spacing  used  in  all  of  these  calculations  except  the 
vertically  averaged  cloud  properties  (which  are  independent  of  the  vertical  direction)  was 
Az  =  2  m.  And  finally,  a  constant  time  step  of  At  =  2  s  was  used  throughout  the  simulations. 
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A  qualitative  comparison  of  the  simulated  dispersing  clouds  at  comparable  times  for  each  of 
the  two  tests  is  given  in  Fig.  (2).  In  this  figure,  the  dispersing  clouds  are  represented  by  “dot- 
plots”  where  each  dot  is  the  location  of  an  individual  marker  particle  projected  onto  the  plane  of 
interest.  In  Fig.  (2a)  the  cloud  is  viewed  from  above  looking  down  onto  the  horizontal  plane  and 
in  Fig.  (2b)  the  cloud  is  viewed  from  the  side  with  the  wind  blowing  from  left  to  right  in  both 
views. 

While  the  gravity  flow  rates  are  similar  in  both  simulations  (on  the  order  of  .5  m/s),  the  effect 
on  cloud  shape  is  more  pronounced  in  the  low  wind  speed,  more  stable  test  of  Burro  8.  When 
viewed  from  above  [Fig.  (2a)],  the  simulated  Burro  8  cloud  is  seen  to  be  rounder  or  “pancake” 
shaped,  while  the  Burro  9  cloud  has  the  more  typical  “cigar”  shape.  In  addition,  the  simulated 
Burro  8  cloud  is  wider  [Fig.  (2a)]  and  lower  [Fig.  (2b)]  than  the  Burro  9  cloud.  Under  the  low 
wind  speed  conditions  of  Burro  8,  gravity  flow  has  more  time  to  widen  the  cloud  before  it  travels 
downwind  and  dissipates  to  where  the  density  is  essentially  equal  to  the  ambient  value. 
Consequently,  gravity  flow  effects  on  cloud  shape  are  greater  in  Burro  8  than  in  the  higher  wind 
speed  case  of  Burro  9.  Similarly,  the  dense-gas  effect  on  cloud  height  is  also  greater  in  the 
simulated  Burro  8  test,  resulting  in  a  lower  cloud  in  the  Burro  8  simulation  in  comparison  to  the 
Burro  9  simulation. 
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Figure  2.  Dot-plot  projections  of  the  marker  particle  locations  onto  the  (A)  horizontal  and  (B)  vertical-inwind  planes  with  the  wind  direction 
pointing  from  left  to  nght  in  both  views.  Plots  (Al)  and  (Bl)  are  top  and  side  views  of  the  Burro  8  simulation  at  4  min.  Plots  (A2)  and  (B2)  are 
top  and  side  views  of  the  Burro  9  simulation  at  1.5  min. 
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Table  2.  A  model-data  comparison  of  peak  concentration  expressed  as  the  volume  fraction. 


Dist.  (m) 

Burro  8 

Calc. 

Obs. 

57 

.63 

.58 

140 

.34 

.25+ 

400 

.11 

.08 

800 

.038 

.024 

Dist.  (m) 

Burro  9 

Calc. 

Obs. 

57 

.60 

140 

.28 

.10 

400 

.074 

.035 

800 

.026 

.012 

These  qualitative  features  of  the  simulated  clouds  agree  with  what  was  observed  in  the  actual 
experiments.  A  more  quantitative  comparison  between  the  simulated  cloud  results  and 
experimental  observation  is  given  in  Table  2  where  the  peak  volume  fraction  is  compared  at  the 
four  downwind  distances  where  measurements  were  made.  Close  to  the  source,  the  calculated 
concentration  is  in  very  good  agreement  with  observation.  At  greater  distances  from  the  source, 
the  calculated  concentration  becomes  more  conservative  (tends  to  over-estimate  concentration)  in 
comparison  to  observation.  At  the  largest  downwind  distance,  the  calculated  concentration  is 
about  50%  greater  than  observed  in  Burro  8  and  100%  greater  in  Burro  9. 

In  a  qualitative  comparison  of  the  cloud  width  and  height  (not  shown),  the  simulated  cloud 
width  was  in  fair  agreement  with  observation,  but  the  simulated  cloud  height  was  somewhat 
lower  than  observed.  This  suggests  that  the  over-estimation  of  cloud  concentration  at  larger 
distances  may  be  due  to  an  under-estimation  of  cloud  height  in  this  region. 
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CONCLUSION 


A  dense  gas  version  of  the  ADPIC  Lagrangian  particle,  advection-diffusion  model  has  been 
developed  that  is  capable  of  dispersion  simulations  over  terrain  with  time-  and  spacially-varying 
ambient  winds.  The  new  dense-gas  version  includes  the  physical  processes  associated  with 
dense-gas  dispersion  that  affect  the  ambient  windfield,  turbulent  diffusivity,  and  thermodynamic 
properties  of  the  dispersing  cloud.  These  effects  result  in  a  lower  and  wider  cloud  than  is 
observed  when  a  trace  or  neutral  density  gas  is  released  and  in  a  cloud  that  lends  to  follow  the 
downhill  slope  of  terrain  independent  of  the  wind  direction.  The  new  dense-gas  version  of 
ADPIC  is  applicable  to  ground  level  releases  of  both  ambient-temperature  and  colder-than- 
ambient,  denser-than-air  gases.  The  dispersion  region  treated  by  the  model  extends  from  the 
near-field  where  dense-gas  effects  are  most  pronounced  out  to  die  far-field  where  the  originally 
dense-gas  cloud  disperses  like  a  trace  gas  cloud. 

During  the  course  of  this  effort,  several  numerical  and  scientific  difficulties  arose  which 
needed  to  be  resolved  in  order  to  create  a  reliable  and  robust  dense-gas  version  of  die  code.  Three 
of  the  major  obstacles  were  (1)  the  inability  of  a  fixed  grid  to  resolve  the  vertical  structure  and 
height  of  low  lying  dense-gas  clouds,  (2)  the  inability  of  advection-diffusion  models  to  simulate 
thermal  transport,  and  (3)  the  occurrence  of  spurious  gravity  waves  due  to  the  unevenness  in  the 
instantaneous  marker  particle  distribution.  The  vertical  resolution  problems  were  solved  by  using 
the  vertical  averaging  technique  and  assumed  vertical  profiles  to  calculate  the  thermodynamic 
properties  of  the  dense-gas  cloud  and  then  using  stochastic  methods  to  calculate  the  dynamics  of 
the  marker  particle  displacements.  The  inclusion  of  thermal  transport  was  added  to  the  model  by 
defining  a  time-varying  “thermal  energy  deficit”  property  for  each  marker  particle.  Then,  as  the 
marker  particles  disperse  with  time,  thermal  transport  is  calculated  from  the  average  “thermal 
energy  deficit”  within  each  region  of  the  cloud.  And  finally,  the  problem  of  spurious  gravity 
waves  was  overcome  by  using  spatial  averaging  of  the  gravity  forcing  functions. 

The  results  of  the  model  simulations  conducted  to  date  compare  favorably  with  experimental 
observation  with  regard  to  the  centerline  concentration,  the  general  cloud  shape  and  the  increase 
in  cloud  width  due  to  crosswind  gravity  flow  of  the  dense-gas  cloud.  The  predicted  concentration 
appears  to  become  somewhat  conservative  (over-predict)  with  downwind  distance  (generally 
within  a  factor  of  two)  and  this  result  appears  to  be  connected  to  under-estimating  the  growth  rate 
of  the  dense-gas  cloud  height.  Possible  causes  include  over-estimating  the  damping  effect  of 
stable  density  stratification  on  the  turbulent  diffusivity,  under-estimating  the  turbulence 
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producing  effect  of  surface  heating,  and  over-estimating  the  gravity  flow  rate  in  the  higher 
elevations  of  the  cloud.  Since  the  discrepancies  between  the  model  predictions  and  experiment 
are  within  a  factor  of  two,  additional  comparisons  are  needed  to  determine  whether  the 
differences  are  due  to  systematic  deficiencies  in  the  model  or  the  normal  variation  of  the 
atmosphere. 

And  finally,  a  secondary  goal  of  this  work,  beyond  developing  a  technique  for  simulating 
dense-gas  dispersion  over  complex  terrain,  was  to  structure  these  new  capabilities  in  a  manner 
that  they  could  be  readily  transferred  to  other  codes.  Toward  this  end,  a  modular  approach  was 
used  to  include  the  dense-gas  modifications  within  ADPIC,  making  them  easier  to  transfer  to 
other  Lagrangian  particle  dispersion  models.  In  addition,  all  of  the  dense-gas  modifications  are 
contained  within  the  ADPIC  dispersion  code,  including  those  that  affect  the  ambient  wind  field 
which  is  calculated  in  a  separate  code.  Consequently,  the  new  dense-gas  version  of  ADPIC  is  not 
restricted  to  a  single  wind  flow  model;  rather,  it  can  be  run  with  alternative  wind  flow  models  to 
provide  the  ambient  wind  field.  These  models  include  the  fully-conservative  (mass,  momentum, 
and  energy)  prognostic  wind  flow  models,  as  well  as  the  more  commonly  used  data- 
interpolating,  mass-conserving,  diagnostic  models.  This  versatility  of  the  dense-gas  version  of 
ADPIC  allows  it  to  be  coupled  with  a  diagnostic  wind  flow  model  for  emergency  response 
purposes,  and  then  when  time  constraints  are  reduced,  it  can  be  used  with  a  physically  more 
complete,  yet  computationally  slower,  prognostic  model. 
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